R语言执行单因素协方差分析
前文简介了单因素方差分析（单因素ANOVA）在R语言中的实现方法，本篇继续简介单因素协方差分析（ANCOVA）。当方差分析中存在协变量时，即可称为协方差分析。其中单因素协方差分析是最常见的，在单因素方差分析中引入了协变量。
何谓协变量，举例来说。对怀孕母鼠喂食某种药物，并观察药物处理组和对照组相比，新生小鼠体重是否具有区别。由于各母鼠的怀孕时间有所差别，而这个怀孕时间可能会对小鼠体重产生影响而不可忽略，就可视作协变量处理。再举一例，重测序分析中统计基因的突变频率，由于各基因长度是不一样的，即使碱基是随机突变的，更长的基因也可能会出现更多的突变碱基数量。此时基因的长度就是一个协变量，必须考虑在内。好了，下面接正文。
本文使用的作图数据的网盘链接（提取码z4w4）：
https://pan.baidu.com/s/1J-9GsmoHuQ_CEpxeWyEQsA
数据预处理
示例数据说明
我们首先将示例数据读到R中，并从中挑选部分数据作为演示。
#读入文件soil <- read.table('soil.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
group <- read.table('group.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
soil <- merge(soil, group, by = 'sample')
#以 shannon 指数为例，同时将分组列转换为因子变量
shannon <- soil[ ,c('sample', 'treat', 'shannon', 'days')]
shannon$treat <- factor(shannon$treat)
str(shannon)
head(shannon)
假设存在这么一个研究：
我们使用来源于同一环境中的土壤进行盆栽试验（土壤类型一致），并种植了同种植株（植物类型一致）。我们将盆栽（1植株/1盆栽）分为了3组，分别在土中添加了三种化学物质（a、b、c）；然后等待植物到达花期后，收集每个植株的根际土，并通过16S测序，获得了植物根际细菌群落的Alpha多样性指数，意在探究不同类型的化学物质是如何影响植物根际细菌群落的。
由于各盆栽中植株到达开花期时间并非一致，无法保证所有植株均在同一天取样观察。尽管期间并未相隔很长时间（天），理论上单因素ANOVA就可以满足需求，但我们仍然想要将各盆栽中植株到达花期的时间作为协变量去考虑（如果觉得植物生长时间相隔几天的差异可能会导致其根际菌群产生较大的变异时，尽管实际上可能有些不妥，但作为示例，请忽略这点），尝试使用单因素协方差分析（ANCOVA）。
对应于上述挑选出的测试数据“shannon”：sample，试验样本名称，每一个样本即对应一个盆栽，各盆栽中土壤类型、植物类型完全相同，而且均是1植株/1盆栽；treat，在土壤中添加的三种化学物质（a、b、c），这列作为分组列，需要转换为因子变量类型，各组处理间相互独立；shannon，Alpha多样性指数中的Shannon指数，数值变量；days，各盆栽中植株的开花时间（即生长天数），这里为数值变量。
这里，植物根际菌群的Shannon指数为因变量，植物生长天数（days）为协变量，考虑使用ANCOVA，探究不同类型的化学物质处理下的植物根际细菌群落的Shannon指数是否存在显著不同。
评估检验的假设条件
ANCOVA要求数据服从正态分布，以及各组方差相等，同时还假定回归斜率相同。
正态性检验
同单因素ANOVA的方法，可使用Q-Q图来检验正态性假设。#QQ-plot 检查数据是否符合正态分布library(car)
qqPlot(lm(shannon~treat, data = shannon), simulate = TRUE, main = 'QQ Plot', labels = FALSE)由图可知，我们的数据符合正态分布模型（尽管似乎有个离群点，这时候可以考虑删除这个样本后再继续，本示例演示不再将它删除，直接进入下一步分析）。
bartlett.test(shannon~treat, data = shannon)结果显示，我们的数据各组方差相等。
回归斜率的同质性检验
这里ANCOVA包含植物生长时间×化学物质类型的交互项，需对回归斜率的同质性进行检验，若交互效应显著，则意味着植物生长时间和植物根际菌群的Shannon指数的关系依赖于所添加的化学物质类型。#检验回归斜率的同质性，交互效应不显著则支持斜率相等的假设（即 p 值大于 0.05 说明回归斜率相等）summary(aov(shannon~days*treat, data = shannon))根据aov()公式可知，这实际上是一个双因素方差分析（关于双因素ANOVA，下篇再简介），根据双因素方差分析中的交互项结果来判断回归斜率的同质性。结果显示交互作用不显著，支持了斜率相等的假设。如果假设不成立，可以尝试变换协变量或因变量，或使用能对每个斜率独立解释的模型，或使用无需回归斜率相等的非参数ANCOVA方法（如sm包sm.ancova()）。
单因素协方差分析（ANCOVA）
单因素ANCOVA
上述检验均已通过，接下来进行单因素ANCOVA。
如果不满足上述前提假设，一是可以考虑转化数据（当然，我们需要确保转换后的数据能够被合理解释，否则将无意义），二是可以考虑使用非参数的检验方法。上述提及了一个对应单因素协方差分析的非参数替代方法，sm包sm.ancova()。
对于带协变量的项，以单因素协方差为例，aov()函数书写为aov(y~x+A)的样式，其中x为协变量，A为因子变量，注意需要将协变量写在因子前面，如上所示，协变量植物生长时间（days）在前，因子化学物质类型（treat）在后，顺序不可颠倒。
#满足假设，单因素协方差分析，详情使用?aov查看帮助fit <- aov(shannon~days+treat, data = shannon)
summary(fit)
ANCOVA结果表明：（1）植物生长时间相隔几天的差异并未导致其根际菌群产生较大的变异（p值不显著）；（2）控制植物生长时间，不同类型的化学物质处理下的植物根际细菌群落的Shannon指数存在显著不同（p值显著）。
对于各组均值的获得方式。
#查看各组均值及标准差aggregate(shannon$shannon, by = list(shannon$treat), FUN = mean)
aggregate(shannon$shannon, by = list(shannon$treat), FUN = sd)
#由于使用了协变量，若想获取去除协变量效应后的组均值（调整的组均值）
library(effects)
effect('treat', fit)
因变量、协变量和因子之间的关系
#详情使用 ?ancova 查看帮助
library(HH)
ancova(shannon~days+treat, data = shannon)
ancova(shannon~days*treat, data = shannon)
如上示例，通过“ancova(shannon~days+treat, data = shannon)”，我们又执行了一次ANCOVA，结果屏幕输出，和上文结果一致。同时通过关系图可知，3条回归线相互平行，只是截距项不同，b组截距项最大，c组截距项最小；回归线拟合效果并不理想，也对应了我们先前的结论，在“数天”这么一个短期时间内，植物根际菌群的Shannon指数没有发生剧烈的改变。
我们还通过“ancova(shannon~days*treat, data = shannon)”，执行了一次双因素ANOVA。不过在这里我们并不是想做个双因素ANOVA分析，而是在更改了函数公式后，意在可视化允许斜率和截距项依据组别而发生改变，从而体现那些违背回归斜率同质性的实例。（上文已知，回归斜率的同质性检验是通过双因素ANOVA中的交互作用判断的，本示例中的回归斜率的同质性检验已通过，大家可以尝试自行寻找一例回归斜率不相等的数据，然后使用ancova()作图查看其与回归斜率相等的数据的差异）
友情链接